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SUMMARY 

Numerical flux functions for solving the Euler equations using exact and/or approximate 
solutions of the Riemann problem of gasdynamics are discussed. Under certain restric- 
tive conditions, schemes using these flux functions produce systems of equations which can 
exhibit a single degree of freedom. In some instances, the solutions represented by this 
degree of freedom are unstable to perturbations. This local instability can seriously degrade 
the temporal convergence of numerical schemes. This point is demonstrated by numerical 
example. 

1. INTRODUCTION 

A major portion of these notes should be viewed as a review of some known (but perhaps 
not well known) technical details of flux functions based on approximate and/or exact 
solutions of the Riemann problem (see van Leer (1984), Roe (1981), Osher and Solomon 
(1981)). The remainder of these notes describes some anomalous behavior resulting from 
these constructions. The most remarkable result appears to be that these anomalies can 
even be found in flux functions which require the exact solution of the Riemann problem. 

2. PRELIMINARIES 

We consider schemes for solving the one-dimensional inviscid flow equations 


u, + f(u), = 0 


( 2 . 1 ) 


where u represents the vector of conserved variables for mass, momentum, and energy. 
Numerical schemes in semi-discrete form will be assumed to take the following form for the 
control volume spanning *j_i < x < Xj+i 


d\ii hi-i- 1 ~ h _ i 

J _l_ jZ i j % _ q 


dt 


Ax 


( 2 . 2 ) 


In this equation, uy is the discrete approximation to u(xj) and h J± ± are the numerical flux 
functions evaluated at the cell-volume interfaces. 

In the present analysis we consider several numerical flux functions for solving equation 

2 . 2 : 


Godunov/Riemann Flux 

h(u L , U R ) = f(u), U = u Riemann{x/t = 0; U L , U R ) (2.3a) 


S-C-S (Shock-Contact-Shockl Riemann Flux 

h(u L ,u R ) = f(u), u = u a . c . s (x/t = 0; u L , u fl ) (2.3b) 


Roe Flux 

= i(f H + f L ) - lsign(^) RM (f a - f 1 ) 


(2.3c) 



State Fluxes 


h(u L ,u«) = f(v), v = i(v R + v L ) - ^sign {M-' A Roe M){v R - v L ) (2.3d) 

In these equations, A is the usual flux Jacobian matrix, A = df/d\\. The S-C-S Riemann 
flux is constructed from elementary wave-transition operators corresponding to compression 
shocks, entropy-violating expansion shocks, and contact discontinuities. For the state flux, 
the matrix M is required to satisfy the mean value construction, Au = M(u R , u L ) Av, and 
pointwise v = M -1 u. This defines a multitude of possible state vectors. For simplicity, we 
will consider only three such constructions for this flux wherein v corresponds to the vector 
of conserved variables [p,pu, pe], vector of primitive variables [/>, u,p], and Roe’s parameter 
vector y/p [l,u, IT]. 

In the next section the underlying shock structure for these fluxes is reviewed. This is 
an important aspect of the analysis because this information allows us to clearly identify 
the underlying degree of freedom (d.o.f.) in the discrete equations when discontinuities are 
present. By using this information, local stability and convergence characteristics can be 
explored. 

3. REVIEW OF NUMERICAL STATIONARY SHOCK STRUCTURES 

In this discussion, the dependent variables are assumed to be piecewise constant and the 
numerical flux is obtained by solving an exact or approximate Riemann problem centered 
at each cell face. Rather than consider arbitrary data, left and right state data, u L and u R , 
are chosen to satisfy the stationary Rankine-Hugoniot jump relations, f L — f R = 0 subject 
to an intermediate state, u M , as shown in figure 1. 
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Figure 1 . Shock transition profile. 

Stationary resolution implies constancy of numerical fluxes. 

h L = h LM = h MR = h R 

For numerical fluxes constructed of the form, h = f(u), this also implies that 

f(u L ) = f(5 LM ) = f(u Mfl ) = f(u H ) 

which requires the simultaneous satisfaction of the stationary jump conditions ( [f]“ b = 
f(u 6 ) - f(u“)) 


[f] 


L-LM 


= 0, [f] 


LM-MR 


= 0, [f] 


MR-R 


(3.1) 


2 




Figure 2. Shock Hugoniot curves. 


This occurs only if u LM and u MR take on either the right or left state values. To see this, 
consider the implications in the p — r plane as shown in figure 2 (r = 1/p). In this figure, 
the Hugoniot curves, H(p,T;p L ,t l ) and H(p,T-,p R ,T R ), depict possible shock solutions 
(moving or stationary) connecting to the left and right states, respectively. The initial 
data is prescribed such that H(p,T-,p L ,r L ) fl H(p,r-,p R ,r R ) at the left and right states. 
Equation 3.1 requires that u LM connect to both left and right states by shock solutions. 
This is satisfied only at the intersection points of the two Hugoniot curves where u LM 
takes on either the left or right state values. To determine which is correct, we look at the 
functional dependency of the flux state, u LM — u(u L , u M ). Because we are interested in 
cases such that u M u R , we conclude that u LM = u L . A similar argument holds for u MR . 
For the class of flux functions of the form, h = f(u), the following conditions are necessary 
for cell equilibrium when u R and u L connect by a single stationary discontinuity 


5 LM 


= u 


~M R _ U -R 


(3.2) 


This restricts the possible values of u M . We will use equation 3.2 to show that the 
numerical scheme, equation 2.2, with Godunov flux function allows resolution of stationary 
discontinuities with one transition point. The parallel theory for flux functions, equations 
2.3b-d, will then be discussed. 


Godunov and S-C-S Flux Subcell Wave Structure 

The requirement that conditions specified in equations 3.2 are satisfied allows one to 
obtain qualitative and quantitative information as to the particular waves present in each 
cell-face Riemann problem as illustrated in figure 3. To gain some insight, we might begin by 
considering the right cell face and its Riemann problem in the limit as u M approaches u R . 
In this case, the (2)-(3) waves travel to the right and the (l)-wave travels to the left with 
speeds approaching their characteristic limits. Before the characteristic limit is reached, the 
waves would have finite amplitude. Yet the Riemann problem at the right cell face must 
still produce a state for x/t — 0 equal to the right state value, u fl . This can happen only if 
the amplitudes of the (2)-(3) waves are identically zero. 

Now relax the condition that u M and u R be close together. This wave structure ((2)- 
(3) waves of zero amplitude) remains valid. If the (l)-wave is a shock, then the jump 
conditions are satisfied, [f] M_ = a [u] ~ and valid states, u M , lie on the Hugoniot 

curve, H(p,r-,p R , t r ) connecting u M to u fl . Adding the condition that the shock propagate 
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Figure 3. Cell wave schematic for shock transition. 


into the cell (a < 0) we have constraints on p, p, and u 



p L <p M < p R 

(3.3a) 


p L <p M < p R 

(3.3b) 


u L > u M > u R 

(3.3c) 

The conditions p M < p R and p M < p R ensure that the (l)-wave is a shock and p L < 
p M and p L < p M ensure that the shock-propagation speed is negative. The bounds on 
velocity follow from the “mechanical” shock conditions (Courant and Friedrichs, 1948), 
G(p,T',p R ,T R ) = (u — u R ) 2 . If the (l)-wave is a rarefaction, then ( u — c) M < (u — c) R < 0 
(the (l)-wave always travels into the cell) and we obtain limits on p, p, and u 


L s' R ^ M 

P < P <P 

(3.4a) 


p L <p R < p M 

(3.4b) 


u L >u R > u M 

(3.4c) 


As yet, we have not guaranteed that conditions, equation 3.2, on the left cell face are 
satisfied. Only by considering the left cell-face Riemann problem will we find that not all 
of these states are possible. 

The Riemann problem structure for the left cell face is more complicated. Examining the 
limit as u M approaches u L , general Riemann problem solutions are possible with all three 
waves of nonzero amplitude. Conceptually this means that global information concerning 
the Riemann problem solution is needed. Recall that the solution of the Riemann problem 
is a composition of elementary wave-transition operators corresponding to shocks, contacts, 
and rarefactions 

U R = T (3) T (2) r (l) U L (3.5) 

Each transition operator can be constructed as a function of a scalar parameter x. Ob- 
taining a solution of equation 3.5 implies finding the three real valued scalars, X\,X 2 , and 
X 3 . Fortunately, we need determine only xi to find all additional constraints on u^. These 


4 




additional constraints will guarantee positive propagation speed of the (l)-wave and, re- 
gardless of the particular (2)-(3) waves present, no additional constraints on u M can be 
imposed. The first step is to determine if the (l)-wave is a shock or rarefaction wave. In 
solving the Riemann problem, equation 3.5, Smoller (1982) reduces the system to solving 
a single nonlinear equation in X\. Once xi is known, *2 and X3 are easily calculated. In 
his notation we have that the (l)-wave component of the Riemann problem solution is a 
rarefaction wave if and only if Xi > 0 is a solution of 


4>(x 1 ) = h(x 1 ) + \j — h(x x + log(R)) = C 


with 


h{x 1 ) = < 


( *(» - 


e~ TXl ), X!>0 


L -e 11 x, < 0 

7-1 (l+/ 3 e-*l)l /2 ’ - u 


P M P M 

A = -T' B =~l' G = 

p L p L 


u M — u L 


T = 


7-1 

2 7 


and is a shock otherwise. Independent of which (l)-wave is present on the right cell face, we 
have, from equations 3.3 and 3.4, A >1,5 > 1, and C < 0. A simple analysis reveals that 
under these conditions the only solutions of 4>{x\) = C are for xi < 0 and the left cell-face 
(l)-wave must be a shock wave. The shock speed in this case is given by 


cr — u L — c L 


(/?-!)* 

( 3 - z 


1/2 


with /? = 7±1 and z — . Positive propagation speed is achieved whenever 


= - log 


'((3 + i)(u L /c L y 


< X! < 0 


(3.6) 


We can readily compute the limit cases of <p(xi) = C 

0, 


x\ = < 


U M = U L 


it i 


U M = U R 


The function <f>(x 1 ) is a monotonically increasing function of its argument and C varies 
monotonically with u M . From this we can conclude that x m ^ n < x\ < 0 only for u L < 
u M < u R . This excludes the right cell-face rarefaction solution. Thus, the only solution 
remaining connects u M to u R on the Hugoniot curve, H(p,r;p R ,r R ), and connects u L to 
in general by a full three-wave solution of which the (l)-wave must be a shock. This 
also guarantees bounds on the variables p M , p M , and u M 

p < p < p 

p L < p M < p R (Godunov/Riemann flux) (3.7) 
u R > u M > u R 
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Without actually knowing the amplitude and/or nature of the (2)-(3) waves on the left 
cell face, we can consider the picture “complete”, i.e., the (2)-(3) waves can add no further 
constraints on . The remainder of this section will be devoted to a similar analysis of the 
flux functions (2.3b-d). Since the resultant bounds will be identical to that of the Godunov 
flux, the reader may wish to proceed directly to section 4 without loss of continuity. 

The analysis for the S-C-S Riemann flux follows a similar set of arguments. The right 
cell-face wave structure connects u M to u R by a single (l)-wave compression or expansion 
shock solution. For the (l)-wave compression shock we have the conventional limits 


P L <P M < P R 

(3.7a) 

P L <P M < P R 

(3.7b) 

u L > u M > u R 

and for a (l)-wave expansion shock we have 

(3.7c) 

P L <P R < P M 

(3.8a) 

P L ~ P R — P M 

(3.8b) 

U L > U R > V.M 

(3.8c) 


As with the Godunov flux, the left cell face also requires knowledge of the global Riemann 
problem solution. Following a development similar to the conventional Riemann problem, a 
global solution to this approximate Riemann problem can be constructed. In this case the 
(l)-wave is an expansion shock if x\ > 0 is a solution of 

#ci) = M*i) + +i°g(-ff)) = c 


with 


h(xi) 


2y/¥ 




7 — 1 (1 + /Je - ® 1 ) 1 / 2 

and is a compression shock otherwise. In this equation we have defined the parameters 


A = 


n L ’ 


B 


c - 


u M - u L 


T = 


7~1 

27 


From equations 3.7 and 3.8, we have, A > 1, B > 1, and C < 0, and all solutions of 
<j>(xi ) ■= C are such that x\ < 0. As with the Godunov flux, the (l)-wave is a compression 
shock. Following the same analysis, it follows that positive propagation of the (l)-wave only 
occurs under the following conditions 


M ^ R 

P < P 
P <P 


6 



u M > u R 

From this constraint the only remaining valid solution connects u M to u R by a single moving 
compression shock and connects to by a three-wave solution of which the (l)-wave 
is a compression shock. This wave structure places bounds on p M , p M , and u M identical to 
the Godunov flux 

P L <P M < P R 

p L < p M < p R (S-C-S Riemann flux) (3.9) 
u L >u M > u R 


Roe and State Flux Subcell Wave Structure 

The analysis of the Roe and state fluxes is somewhat simpler. These fluxes also have 
valid intermediate states that connect to the right state by a single moving shock wave. In 
this discussion we will demonstrate only that this wave configuration does indeed satisfy 
the conditions for cell equilibrium. From this we will obtain bounds on the u M state. 
Note that from the shock-jump relations and the Roe identities we have (A Roe = XAJf -1 , 
A = diag(Ai,A 2 ,A 3 ) ) 


<7 [u] = [f] = A Roe [u] 

aX- 1 [u] =AX~ 1 [u] (3.10) 

a [w] =A [w], [w] = X~ l [u] 

In this form the jump relations are satisfied if and only if for some j 

= 0 i = 1, 2, 3 i ^ j 

x . . (A 11 ) 

a = Xi i — j 

From this we have that sign(<r) [w] = sign(A) [w] and the following useful relations: 


sign(cr) [u] — sign(A) j R oe [u] 

sign ( cr) [v] =sign (M* 1 AM) Roe [v] (3.12) 

sign ( <r) [f] =sign(>l)ii oe [f] 

When u M connects to u R by a single moving shock, we have, from equation 3.11, that 
(u — c)“* = <7 and from the condition that cr < 0 we can conclude 


P L <P M 

(3.13a) 

% 

VI 

(3.13b) 

Al 

(3.13c) 


The flux formulas (equations 2.3c-d) can now be shown to satisfy the right cell-face re- 
quirements (equation 3.2) subject to the constraints (equations 3.13a-c). In this case, 
sign(cr AfR ) = —1, and from equations 3.12 we have 
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Roe flux (right cell face): 


h MH = 1 (f M + fR) _ l sign (^) floe (fH _ fM) = 1 (f M + f fl) + I(f« - f M ) 
=f R 


I(v M -f v R ) + ^(v fl — v M ) 

Z Z 

We now consider the left cell face. For Roe’s flux we require 

i L = h LM = + i L ) - ^sign(A) fioe (f A/ - f L ) (3.14) 

After some rearrangement we obtain 

[f] L - M = sign^)*,. [f] 1 -" 

or in terms of the wave amplitudes 

[w] L-M = sign(A) [w] L-M (3.15) 

This equation would be satisfied if u M = u R and u L and u M connect by a stationary shock 
(with sign(O) = 1) or whenever [Aj , ^ 2 , A 3 ] > 0 . The former situation is actually a limiting 
case and both are encompassed by the condition [Ai,A 2 ,A 3 ] iM > 0 . This amounts to 
showing that 

(« - < 0*2 > 0 ( 3 - 16 ) 

subject to constraint that u M lies on H(p,T;p R ,r R ) = 0 and satisfies the constraints listed 
in equation 3.13 . The analysis for the state fluxes would be identical and will not be repeated 
here. Although we do not prove this here, we find that these conditions are satisfied only 
when 

p M <p R 

p M < p R 

u M > u R 

Combining these conditions with equation 3.13 we obtain the now familiar conditions 

p < p < p 

p L < p M < p R (Roe and state fluxes) (3.17) 

u L > u M > u R 


State fluxes (right cell face): 
h MR =f(u MR ) = f(u(v MH )) 

\ MR =^(v M + v R ) — ^sign(Af AM)R oe (v R — v M ) = 

Z Z 
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Under these conditions, we have that sign(A) = I and the requirements for cell equilibrium 
are satisfied 

Roe Flux (left cell face): 

h LM = 1 (f L + f M) _ l sign(j4 ) Hoe(f M _ {L) = 1 (f L + f M } _ 1 (f M _ f L } 

=f L 


State Fluxes (left cell face): 

h LM =f(u iM ) = f(u(v LM )) 

v LM =J(v L + v M ) - isign (M - 1 AM) Roe (v M - v L ) - i(v L + v M ) - i(v M - v L ) 
Z 2* L L 


4. STATIONARY SHOCK WAVE CHARACTERISTICS 

We have shown that all fluxes considered require that the primitive variables be bounded 
between the left and right states 

p L <p M < p R 

P L <P M <P R (4.0) 



Figure 4. Transition point momentum deviation for M=10 shock. 


Note that interval boundedness in primitive variables does not guarantee interval bound- 
edness in other variables. Two obvious examples are momentum and stagnation enthalpy. 
Across a stationary shock these quantities are constant, yet values at the transition point 
may deviate significantly 


p L u L = p R u R / p M u M 


(4.1a) 
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H l = H r ± H m (4.1b) 

In figure 4 we demonstrate this phenomenon for data corresponding to a Mach 10 stationary 
shock. In this case we find that, depending on the particular location of the state, u M , 
we can have up to a 40% spike in momentum. This can adversely affect some numerical 
schemes which interpolate/extrapolate state variables and/or other quantities. In part 2 of 
these notes (T.J. Barth, in preparation), we examine this effect in detail in the context of 
propagating shock waves. 

The wave structure of these fluxes clearly identifies a single d.o.f. present when stationary 
discontinuities are present. For example, if p M is chosen as the free parameter, then any 
value can be chosen subject to P L <P M < p R from which the entire state, u M , is fixed such 
that cell equilibrium is satisfied. Pressure is obtained from H(p,r]p R ,t r ) = 0 and velocity 
from G(p,r;p R ,T R ) = ( u R — u ) 2 . This d.o.f. implies that the residual function 

r(u M ; u L , u R ) = h MH — h LM (4-2) 

has a singular Jacobian matrix, i.e., det(^^-) = 0. To see this, construct a vector in state 
space (p,u,p) tangent to the curves, H(p,r',p R , t r ) = 0 and G(p,T',p R ,r R ) = ( u R — u ) 2 . 
Then construct its image, v, in ( p,pu,pe ) space. From the definition of a Frechet derivative 
it can be seen that this vector is an eigenvector associated with a zero eigenvalue 

lim -(r(u M + ev) — r(u M )) = v = 0 (4.4) 

o e v ” du M 

To understand the consequences of this zero eigenvalue, one might consider solving a system 
with one interior point, u M , and boundary data, u L , u R 


ft „M 

+ r(u M ;u L ,u H ) = 0 (4.5) 

Near a stationary solution, r(u*;u £ ',u B ) = 0, the system can be locally linearized and the 
solution error, e = u M — u*, behaves according to the equation 



dv 


(u*) e = 0 


(4.6) 


Let x* and A; denote the ith eigenvector and eigenvalue of dr/du M at u*. The solution of 
equation 4.6 is given by 

3 

e (*) = XI c i x * e_Ait 

i= 1 

where Cj are determined from e(0). Stability of this equation requires that A; > 0. The 
presence of a zero eigenvalue implies that if the solution is perturbed from a fixed point, u*, 
the solution need not return to the same fixed point, but can choose another according to 
the equation 

3 

e(f) = cjXj + ^CiXie -A,t 
»=2 
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For the locally linearized system, we have, from equation 4.4, that xj = v and the solution 
is free to move along a local tangent to the shock curves. Note that the singularity of 
the system does not affect the rate of convergence but merely reflects a d.o.f. Note that 
this singularity does affect algorithms such as Newton’s method for solving r(u) = 0 which 
require invertibility of the Jacobian matrix. Also note that if we had considered a system 
consisting of severed computational cells in space with a stationary shock somewhere in the 
interior, the resulting Jacobian matrix would still be exactly singular at the fixed point. 
The following matrix represents this situation (Tj = dr/du M , |A| = A + — A ~ ): 


A+ L 

-a+ l a+ l 

-A +l T\ t 2 

-Tj T 3 A~ r 

-a +r \a\ r a~ r 

- a+ r \a\ r A- r 

-A +r \a\ r 


This matrix has an obvious partitioning which is block lower triangular and to determine 
whether this matrix is singular we need only examine the following submatrix: 

T 2 

t 3 A- r 
a +r \a\ r A- r 

-a +r \a\ r A- r 

-A +r |A| fi 

Expanding the determinant of this matrix about the first column, we can easily see that 
this matrix is singular whenever T\ = dr/du M is singular. 

5. SHOCK TRANSITION POINT STABILITY 

We have shown which shock structures are possible for the flux functions, equations 2.3a- 
d; we now check the stability of these solutions to perturbations. Stability requires that the 
eigenvalues of dr/du evaluated at the stationary solution are nonnegative. All the fluxes 
considered have the same one-parameter family of possible solutions, i.e., connect u M to u H 
by a single moving shock wave such that the primitive variables are bounded between the left 
and right states. In the following paragraphs we consider a family of stationary shocks with 
preshock Mach numbers ranging from 1.5 to 10.0. We then consider varying a parameter, 
t = 1/p, between its limits. The eigenvalue spectrum of the Jacobian matrix, dr/du, at 
each state is then calculated. In figures 5,6, and 7 this analysis is carried out for systems 



T i 
-T, 
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with 1, 3, and 9 cells, respectively, with fixed boundary data, u L and u R . We first verify 
that one eigenvalue is zero, then we examine the remaining eigenvalues. In particular the 
minimum eigenvalue (real part) for each case is plotted. As can be seen, the results for the 
3 and 9-cell systems are virtually identical, which indicates that all the minimum eigenvalue 
information is generated in the immediate vicinity of the transition point. Calculations with 
many more cells have been carried out with identical results. The results for the single-cell 
calculation differ significantly for all flux functions. This indicates that the interaction with 
neighboring cells is apparently an important factor in the local stability. The results in 
figures 5a-f are presented only to show the need for multiple cells in the calculations. The 
results in figures 7a-f demonstrate the virtual insensitivity to the number of cells (when 
greater than or equal to 3). We therefore can refer all discussion to figures 6a-f. These plots 
indicate clearly that for sufficiently strong shocks all flux formulae have regions in which 
shock transition points are unstable to perturbations. (These occur whenever \ m in < 0.) 
From these graphs we make the following observation: 

For perfect gases (7 = 7/5), all fluxes considered may have transition states u m which are 
unstable to perturbations when the preshock Mach number is greater than about 6. 

It is interesting to note that although the different curves take on different shapes de- 
pending on the different flux formula and shock strength, all fluxes cross the zero axis at 
roughly the same shock strength and values of r. Since the only common feature of these 
fluxes is the representation of the wave structure on the right cell face, we conjecture that 
this may be the ultimate source of these local instabilities. Intuitively this wave structure 
appears to be the most precarious and the requirement that the (2)-(3) waves have zero 
amplitude becomes severe as u M approaches u £ . When the cell is not in perfect equilib- 
rium, these waves have nonzero amplitude and send wave information to neighboring cells 
necessary to drive the system toward equilibrium. It appears this process becomes unstable 
for sufficiently strong shocks. 

This phenomenon has important consequences for unsteady calculations involving propa- 
gating shocks as well as steady-state calculations. The implications for transient calculations 
will be presented in part 2 (T.J. Barth, in preparation). The moving shock is related to 
the stationary shock by a Galilean transformation. In the ideal situation in the limit as the 
shock speed approaches zero (cr — ► 0), the transition point approaches a trajectory which 
moves along the curve H(p,T-,p R ,T R ) = 0 from n R to u L as the shock propagates through a 
given cell. In the state space and time domain only one trajectory is possible that does not 
give rise to spuriously generated waves in space. A necessary condition appears to be that 
all solutions on the curve H(p,T-,p R ,T R ) = 0 for ( p,u,p ) between the pre and postshock 
state must be at least stable in the stationary case. Even this is by no means enough to 
guarantee that spurious waves are not generated. For steady-state calculations, the results 
presented in figure 6 indicate that serious problems may arise. Because the local stability 
characteristics are continuously related to the location of the transition-point, convergence 
characteristics are strong functions of the transition point location on the shock Hugoniot 
curve. If the transition-point trajectory is sufficiently close to a region on the Hugoniot 
curve which exhibits local stability, then the state-space trajectory simply attracts to the 
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3n 



(a) Godunov/Riemann flux. 



(c) Roe flux. 



(b) S-C-S Riemann flux. 



(d) State flux, conserved variables. 



(e) State flux, primitive variables. 



(f) State flux, parameter vector. 


Figure 5. Shock wave profile stability - 1 cell (7 = 7/5). 
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T M / X L 


t m/ t l 


(a) Godunov/ Riemann flux. 


(b) S-C-S Riemann flux. 



(c) Roe flux. 


(d) State flux, conserved variables. 



t“/t l t m /x l 


(e) State flux, primitive variables. (f) State flux, parameter vector. 

Figure 6 . Shock wave profile stability - 3 cell (7 = 7/5). 
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Figure 8. Roe flux convergence histories under two initial conditions. 

Hugoniot curve. If the transition-point trajectory is sufficiently close to the region of the 
Hugoniot curve that exhibits local instability, then the state space trajectory repels from the 
Hugoniot curve and follows some path, which eventually attracts in a stable region. This 
process may or may not be much slower than the former process. If the transition-point 
trajectory is sufficiently close to a region of the Hugoniot curve where |A m j n | < e << 1, 
then convergence can be extremely slow. For example, if Euler explicit time advancement is 
used for solving equation 2.2, then the solution converges at a rate proportional to powers 
of (1 — e). This is demonstrated in figure 8 for a Mach 8 stationary shock computed using 
Roe’s flux and Euler explicit time advancement on a 19-cell domain. The curves in this 
figure graph the £2 norm of the residual as a function of iteration number. After consulting 
figure 6c, two sets of initial conditions where chosen which nearly satisfy the stationary 
residual equations. One set of conditions was chosen with the transition point initialized 
near a region on the Hugoniot curve where A m ; n ~ 1. In this case convergence is rapid with 
the residual dropping 10 decades in about 150 time steps. In the second case the transition 
point is initialized near a region on the Hugoniot curve where X m in vanishes. Using the same 
time step as before, over 2000 time steps were required for the same reduction in residual 
norm. By more careful choice of initial conditions, we have performed calculations requiring 
over 50000 time steps for the same reduction. 

The critical Mach number behavior shows a strong dependence on the ratio of specific 
heats, 7. In figure 9 this dependence is plotted for several values of 7 using the analysis given 
earlier. Note that for 7 — + 1, the critical Mach number approaches 3 and for monatomic 
gases (7 = 5/3), the critical Mach number is apparently infinite. 

6. DISCUSSION 

The preceding analysis explains mechanisms which give rise to extremely slow temporal 
convergence. All fluxes considered indicate that for stationary shocks (with preshock Mach 
numbers greater than about 6) extremely slow convergence may result. Although we have 
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Figure 9. Critical Mach number dependence on gamma. 

delayed discussion of slowly propagating shock waves to part 2 of these notes, we note that 
the stationary properties discussed here can apparently be directly related to anomalies seen 
in unsteady calculations with propagating shock waves. One consequence discussed in part 
2 is the “denting” and “bulging” of slowly propagating shock fronts. 

The results presented indicate that renewed interest in flux functions for strong shock 
waves is needed. 
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